function F = cal_chg_ge(eqm, firm, eqm_cf, firm_cf)
%% Counterfactual changes by ex-ante quintile 

    % Indicators for ex-ante quintile
    Qindicator = eqm.Qindicator; %Nf by 5 matrix
    
    
    % Useful functions
    my_fcn1 = @(x, density) sum(kron(ones(1,5),x.*density).*Qindicator,1)./sum(kron(ones(1,5),density).*Qindicator,1)*100;
    my_fcn2 = @(x, density) sum(kron(ones(1,5),x.*density).*kron(ones(2,1),Qindicator),1)./sum(kron(ones(1,5),density).*kron(ones(2,1),Qindicator),1)*100; 
    my_fcn3 = @(x, density) sum(kron(ones(1,5),x.*density).*kron(ones(3,1),Qindicator),1)./sum(kron(ones(1,5),density).*kron(ones(3,1),Qindicator),1)*100; 

    % Quality
    logdiff = log(firm_cf.qlevel) - log(firm.qlevel);  
    ge_cf.eqt_qlevel_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_qlevel_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_qlevel_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_qlevel_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_qlevel_all           = my_fcn1(logdiff, eqm_cf.density_all);

    % Productivity
    logdiff = firm_cf.lnz - firm.lnz;  
    ge_cf.eqt_z_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_z_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_z_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_z_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_z_all           = my_fcn1(logdiff, eqm_cf.density_all);   
    
    % Wage
    logdiff = log(firm_cf.wage) - log(firm.wage);   
    ge_cf.eqt_wage_switcher    = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_wage_exp_ctn     = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_wage_exp         = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);  
    ge_cf.eqt_wage_nonexp_ctn  = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_wage_all         = my_fcn1(logdiff, eqm_cf.density_all);
    
    % Sales
    logdiff_exp_ctn     = firm_cf.lnsales_exp - firm.lnsales_exp;
    logdiff_nonexp_ctn  = firm_cf.lnsales_nonexp - firm.lnsales_nonexp;
    logdiff_switcher    = firm_cf.lnsales_exp - firm.lnsales_nonexp;  
    ge_cf.eqt_sales_switcher   = my_fcn1(logdiff_switcher,    eqm_cf.density_switcher);
    ge_cf.eqt_sales_exp_ctn    = my_fcn1(logdiff_exp_ctn,     eqm_cf.density_exp_ctn);
    ge_cf.eqt_sales_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_sales_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn,  eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_sales_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);
    
    % Number of suppliers
    logdiff_exp_ctn     = log(firm_cf.indegree_exp) - log(firm.indegree_exp);
    logdiff_nonexp_ctn  = log(firm_cf.indegree_nonexp) - log(firm.indegree_nonexp);
    logdiff_switcher    = log(firm_cf.indegree_exp) - log(firm.indegree_nonexp);   
    ge_cf.eqt_indegree_switcher   = my_fcn1(logdiff_switcher,     eqm_cf.density_switcher);
    ge_cf.eqt_indegree_exp_ctn    = my_fcn1(logdiff_exp_ctn,      eqm_cf.density_exp_ctn);
    ge_cf.eqt_indegree_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_indegree_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn,   eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_indegree_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);

    % Number of customers
    logdiff_exp_ctn     = log(firm_cf.outdegree_exp) - log(firm.outdegree_exp);
    logdiff_nonexp_ctn  = log(firm_cf.outdegree_nonexp) - log(firm.outdegree_nonexp);
    logdiff_switcher    = log(firm_cf.outdegree_exp) - log(firm.outdegree_nonexp); 
    ge_cf.eqt_outdegree_switcher   = my_fcn1(logdiff_switcher,    eqm_cf.density_switcher);
    ge_cf.eqt_outdegree_exp_ctn    = my_fcn1(logdiff_exp_ctn,     eqm_cf.density_exp_ctn);
    ge_cf.eqt_outdegree_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_outdegree_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn,  eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_outdegree_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);
        
%{    
    % Number of workers
    logdiff_exp_ctn     = log(firm_cf.employment_exp) - log(firm.employment_exp);
    logdiff_nonexp_ctn  = log(firm_cf.employment_nonexp) - log(firm.employment_nonexp);
    logdiff_switcher    = log(firm_cf.employment_exp) - log(firm.employment_nonexp);  
    ge_cf.eqt_employment_switcher   = my_fcn1(logdiff_switcher,   eqm_cf.density_switcher);
    ge_cf.eqt_employment_exp_ctn    = my_fcn1(logdiff_exp_ctn,    eqm_cf.density_exp_ctn);
    ge_cf.eqt_employment_exp        = my_fcn2([logdiff_exp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_switcher]);
    ge_cf.eqt_employment_nonexp_ctn = my_fcn1(logdiff_nonexp_ctn, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_employment_all        = my_fcn3([logdiff_exp_ctn; logdiff_nonexp_ctn; logdiff_switcher], [eqm_cf.density_exp_ctn; eqm_cf.density_nonexp_ctn; eqm_cf.density_switcher]);
%}    
       
    % Average log wage of suppliers
    logdiff = firm_cf.avg_wgt_lnwageS - firm.avg_wgt_lnwageS;
    ge_cf.eqt_lnwageS_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_lnwageS_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_lnwageS_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_lnwageS_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_lnwageS_all           = my_fcn1(logdiff, eqm_cf.density_all);    
    
%{    
    % Fraction of higher quality suppliers
    logdiff = log(firm_cf.frac_S_higher_q) - log(firm.frac_S_higher_q);   
    ge_cf.eqt_frac_S_higher_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_frac_S_higher_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_frac_S_higher_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_frac_S_higher_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_frac_S_higher_all           = my_fcn1(logdiff, eqm_cf.density_all);    
    
    
    % Fraction of higher quality customers
    logdiff = log(firm_cf.frac_B_higher_q) - log(firm.frac_B_higher_q);   
    ge_cf.eqt_frac_B_higher_switcher      = my_fcn1(logdiff, eqm_cf.density_switcher);
    ge_cf.eqt_frac_B_higher_exp_ctn       = my_fcn1(logdiff, eqm_cf.density_exp_ctn);
    ge_cf.eqt_frac_B_higher_exp           = my_fcn1(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    ge_cf.eqt_frac_B_higher_nonexp_ctn    = my_fcn1(logdiff, eqm_cf.density_nonexp_ctn);
    ge_cf.eqt_frac_B_higher_all           = my_fcn1(logdiff, eqm_cf.density_all);    
%}    

   
    % Export probability
    diff = firm_cf.prob_exp - firm.prob_exp;  
    ge_cf.eqt_exprob_all = my_fcn1(diff, eqm_cf.density_all);      
    
    % Import probability
    diff = firm_cf.prob_imp - firm.prob_imp;  
    ge_cf.eqt_improb_all = my_fcn1(diff, eqm_cf.density_all);      
    
    % Export intensity
    diff = firm_cf.exp_intst_exp - firm.exp_intst_exp;  
    ge_cf.eqt_exintst_exp = my_fcn1(diff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher);
    
    % Import intensity
    diff = firm_cf.imp_intst_imp - firm.imp_intst_imp;  
    ge_cf.eqt_imintst_imp = my_fcn1(diff, eqm_cf.density_imp_ctn+eqm_cf.density_switcher_im);
    
    
    
%% Overall counterfactual changes
       
    % Useful function
    my_avg = @(x, density) sum(x.*density)./sum(density);

    % Quality
    ratio = firm_cf.qlevel./firm.qlevel;
    ge_cf.overall_qlevel_switcher   = (my_avg(ratio, eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_qlevel_exp_ctn    = (my_avg(ratio, eqm_cf.density_exp_ctn)-1)*100;
    ge_cf.overall_qlevel_exp        = (my_avg(ratio, eqm_cf.density_exp_ctn+eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_qlevel_nonexp_ctn = (my_avg(ratio, eqm_cf.density_nonexp_ctn)-1)*100;
    ge_cf.overall_qlevel_all        = (my_avg(ratio, eqm_cf.density_all)-1)*100;
       
    % Wage
    ratio = firm_cf.wage./firm.wage;
    ge_cf.overall_wage_switcher   = (my_avg(ratio, eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_wage_exp_ctn    = (my_avg(ratio, eqm_cf.density_exp_ctn)-1)*100;
    ge_cf.overall_wage_exp        = (my_avg(ratio, eqm_cf.density_exp_ctn+eqm_cf.density_switcher)-1)*100;
    ge_cf.overall_wage_nonexp_ctn = (my_avg(ratio, eqm_cf.density_nonexp_ctn)-1)*100;
    ge_cf.overall_wage_all        = (my_avg(ratio, eqm_cf.density_all)-1)*100;    
    
    % Average log wage of suppliers
    logdiff = firm_cf.avg_wgt_lnwageS - firm.avg_wgt_lnwageS;
    ge_cf.overall_lnwageS_switcher      = my_avg(logdiff, eqm_cf.density_switcher)*100;
    ge_cf.overall_lnwageS_exp_ctn       = my_avg(logdiff, eqm_cf.density_exp_ctn)*100;
    ge_cf.overall_lnwageS_exp           = my_avg(logdiff, eqm_cf.density_exp_ctn+eqm_cf.density_switcher)*100;
    ge_cf.overall_lnwageS_nonexp_ctn    = my_avg(logdiff, eqm_cf.density_nonexp_ctn)*100;
    ge_cf.overall_lnwageS_all           = my_avg(logdiff, eqm_cf.density_all)*100;
    


%% Return

    F = ge_cf;

    
end